Contexte :
Nous somme employé d’Enercoop, jeune entreprise spécialisé dans la fourniture d’énergie renouvelable
But de l’étude :
Prédire la consommation électrique sur 1 an
Crée un modèle de prévision
Pour cela, nous allons utilisé plusieurs modèles, plusieurs méthodes, de prédiction notamment :
Lissage exponentielle
SARIMA
Machine learning
Prophet
Bibliographie :
## Install library
library(tidyverse)
library(lubridate)
library(fabletools)
library(fable)
library(tsibble)
library(feasts)
library(car)
library(dygraphs)
library(broom)
Nous restreignons l’analyse à la région Grand Est avec comme capteur Strasbourg Entzheim.
La période s’étend 2013 et 2019
Source https://www.rte-france.com/eco2mix/telecharger-les-indicateurs
## Ajout de toutes les années
data_conso_grandest <- rbind(RTE2019,RTE2018,RTE2017,RTE2016,RTE2015,RTE2014,RTE2013)
## modification du df
conso_daily_grandest <- data_conso_grandest %>%
filter(Consommation > 0) %>% # suppression des valeurs a 0
mutate(Date_only = ymd(Date)) %>% # conversion en year month day
mutate(Date_only = as_date(Date_only)) %>% # conversion au format date
mutate(Renouvelable = (Eolien+Hydraulique+Solaire+Pompage+Bioénergies)) %>% # Creation du variable energie renouvelable
select(Date_only, Consommation,Renouvelable) %>% # Selection des variables d'interet
group_by(Date_only) %>% # group by daily avec moyenne
summarise_all(funs(mean(.,na.rm=TRUE))) %>%
as_tsibble(index= Date_only) %>% # conversion du df en tsibble
ungroup()
## Creation d'un dataframe month_avg
conso_month_grandest <- conso_daily_grandest %>%
index_by(yearmonth(Date_only)) %>% # group_by(index_by) month
summarise_all(funs(mean(.,na.rm=TRUE)))
DataExplorer::plot_intro(conso_daily_grandest)
questionr::describe(conso_daily_grandest)
## [2556 obs. x 3 variables] tbl_ts tbl_df tbl data.frame
##
## $Date_only:
## Date: 2013-01-01 2013-01-02 2013-01-03 2013-01-04 2013-01-05 2013-01-06 2013-01-07 2013-01-08 2013-01-09 2013-01-10 ...
## min: 2013-01-01 - max: 2019-12-31 - NAs: 0 (0%) - 2556 unique values
##
## $Consommation:
## numeric: 4672.10416666667 5591.9375 6086.375 5956.52083333333 5408.29166666667 5086.875 6042.875 6619.3125 6852.70833333333 6856.6875 ...
## min: 3295.33333333333 - max: 7849.14583333333 - NAs: 0 (0%) - 2534 unique values
##
## $Renouvelable:
## numeric: 2083.75 1396.75 1644.20833333333 1317.85416666667 1225.04166666667 1166.625 1166.29166666667 1133.875 1163.5625 1171.91666666667 ...
## min: 349.270833333333 - max: 4020.5625 - NAs: 0 (0%) - 2522 unique values
conso_daily_grandest %>%
psych::pairs.panels(
method = "pearson", # correlation method
hist.col = "#00AFBB",
density = TRUE, # show density plots
ellipses = TRUE # show correlation ellipses
)
Pas de données manquante ou de ligne incomplète
2556 observations = 7 ans
La consommation au fil du temps a une moyenne stable = saisonnalité additive
temp_strasbourg_2013_2019 <- mutate(temp_strasbourg_2013_2019,DATE = ymd(DATE)) %>% # conversion en year month day
select(-SNWD) %>% # suppression de la colonne snow
as_tsibble(index="DATE") # conversion du df en tsibble
temp_strasbourg_2013_2019 %>% questionr::describe()
## [2556 obs. x 6 variables] tbl_ts tbl_df tbl data.frame
##
## $NAME:
## character: "strasbourg_entzheim" "strasbourg_entzheim" "strasbourg_entzheim" "strasbourg_entzheim" "strasbourg_entzheim" "strasbourg_entzheim" "strasbourg_entzheim" "strasbourg_entzheim" "strasbourg_entzheim" "strasbourg_entzheim" ...
## NAs: 0 (0%) - 1 unique values
##
## $DATE:
## Date: 2013-01-01 2013-01-02 2013-01-03 2013-01-04 2013-01-05 2013-01-06 2013-01-07 2013-01-08 2013-01-09 2013-01-10 ...
## min: 2013-01-01 - max: 2019-12-31 - NAs: 0 (0%) - 2556 unique values
##
## $PRCP:
## numeric: 1 0 0 0.8 0 0 0 0 0.4 0.6 ...
## min: 0 - max: 66.3 - NAs: 0 (0%) - 153 unique values
##
## $TAVG:
## numeric: 6.6 3.6 3.2 7.1 7.2 6.7 4.7 3.4 2.2 2.9 ...
## min: -7.9 - max: 29.6 - NAs: 0 (0%) - 324 unique values
##
## $TMAX:
## numeric: 8.9 7.2 6.1 8.7 7.9 8.4 9.2 3.8 2.5 4.2 ...
## min: -4.9 - max: 38.9 - NAs: 0 (0%) - 396 unique values
##
## $TMIN:
## numeric: 3.8 0.3 -0.8 5.2 6.3 6 0.2 2.2 1.9 1.3 ...
## min: -11.6 - max: 23.8 - NAs: 0 (0%) - 290 unique values
temp_strasbourg_2013_2019 %>% DataExplorer::plot_intro()
temp_strasbourg_2013_2019 %>% psych::pairs.panels(
method = "pearson", # correlation method
hist.col = "#00AFBB",
density = TRUE, # show density plots
ellipses = TRUE # show correlation ellipses
)
Ici aussi, pas de problème particulier dans les données
Source : https://cegibat.grdf.fr/simulateur/calcul-dju
## Imports données (Grand est 67)
dju_month_grandest <- read.csv("data/calcul_DJU_05_05_2021.csv",dec=",")
dju_month_grandest <- dju_month_grandest %>%
mutate (date = ym(DATE)) %>%
relocate(DJU, .after = last_col()) %>%
as_tsibble(index = date)
En utilisation la méthode de calcul Méteo de Cegibat
seuil_temp = 18 # seuil point de changement température
dju_daily_grandest <- temp_strasbourg_2013_2019 %>%
mutate(
dju_calc_chauffage = case_when( ## calcul DJU chauffage
seuil_temp <= ((TMIN + TMAX)/2) ~ 0,
seuil_temp > ((TMIN + TMAX)/2) ~ seuil_temp - ((TMIN + TMAX)/2),
TRUE ~ 0),
dju_calc_clim = case_when( ## Calcul DJU climatisation
seuil_temp >= ((TMIN + TMAX)/2) ~ 0,
seuil_temp < ((TMIN + TMAX)/2) ~ ((TMIN + TMAX)/2) -seuil_temp,
TRUE ~ 0 ),
dju_calc_total = dju_calc_chauffage + dju_calc_clim # calcul DJU total
)
dju_month_grandest_calc <- dju_daily_grandest %>%
index_by(yearmonth(DATE)) %>%
summarise_at(vars(dju_calc_chauffage:dju_calc_total),~sum(.))
autoplot(dju_month_grandest_calc,dju_calc_chauffage,color='gray')+
geom_line(aes(y=dju_month_grandest$DJU), colour = "#D55E00")+
labs(
y = "DJU",
title = "Comparaision calcul vs cegibat"
)
Les graphiques sont proches, nous pouvons valider la méthode de calcul
## merge des dataframmes
data_est_day <- left_join(conso_daily_grandest, dju_daily_grandest, by = c('Date_only'='DATE')) %>%
select(-NAME) %>% # Suppresion variable name
as_tsibble(index = Date_only)
## Creation d'un data month
data_est_month <- data_est_day %>% # group_by
index_by(yearmonth(Date_only)) %>%
summarise(
across(c(TAVG:TMIN),~mean(.)),
across(c(Consommation:Renouvelable,dju_calc_chauffage:dju_calc_total,PRCP),~sum(.))
)
data_est_day %>% questionr::describe()
## [2556 obs. x 10 variables] tbl_ts tbl_df tbl data.frame
##
## $Date_only:
## Date: 2013-01-01 2013-01-02 2013-01-03 2013-01-04 2013-01-05 2013-01-06 2013-01-07 2013-01-08 2013-01-09 2013-01-10 ...
## min: 2013-01-01 - max: 2019-12-31 - NAs: 0 (0%) - 2556 unique values
##
## $Consommation:
## numeric: 4672.10416666667 5591.9375 6086.375 5956.52083333333 5408.29166666667 5086.875 6042.875 6619.3125 6852.70833333333 6856.6875 ...
## min: 3295.33333333333 - max: 7849.14583333333 - NAs: 0 (0%) - 2534 unique values
##
## $Renouvelable:
## numeric: 2083.75 1396.75 1644.20833333333 1317.85416666667 1225.04166666667 1166.625 1166.29166666667 1133.875 1163.5625 1171.91666666667 ...
## min: 349.270833333333 - max: 4020.5625 - NAs: 0 (0%) - 2522 unique values
##
## $PRCP:
## numeric: 1 0 0 0.8 0 0 0 0 0.4 0.6 ...
## min: 0 - max: 66.3 - NAs: 0 (0%) - 153 unique values
##
## $TAVG:
## numeric: 6.6 3.6 3.2 7.1 7.2 6.7 4.7 3.4 2.2 2.9 ...
## min: -7.9 - max: 29.6 - NAs: 0 (0%) - 324 unique values
##
## $TMAX:
## numeric: 8.9 7.2 6.1 8.7 7.9 8.4 9.2 3.8 2.5 4.2 ...
## min: -4.9 - max: 38.9 - NAs: 0 (0%) - 396 unique values
##
## $TMIN:
## numeric: 3.8 0.3 -0.8 5.2 6.3 6 0.2 2.2 1.9 1.3 ...
## min: -11.6 - max: 23.8 - NAs: 0 (0%) - 290 unique values
##
## $dju_calc_chauffage:
## numeric: 11.65 14.25 15.35 11.05 10.9 10.8 13.3 15 15.8 15.25 ...
## min: 0 - max: 26.25 - NAs: 0 (0%) - 548 unique values
##
## $dju_calc_clim:
## numeric: 0 0 0 0 0 0 0 0 0 0 ...
## min: 0 - max: 11.2 - NAs: 0 (0%) - 226 unique values
##
## $dju_calc_total:
## numeric: 11.65 14.25 15.35 11.05 10.9 10.8 13.3 15 15.8 15.25 ...
## min: 0 - max: 26.25 - NAs: 0 (0%) - 603 unique values
data_est_day %>% DataExplorer::plot_intro()
data_est_day %>% psych::pairs.panels(
method = "pearson", # correlation method
hist.col = "#00AFBB",
density = TRUE, # show density plots
ellipses = TRUE # show correlation ellipses
)
autoplot(data_est_day,Consommation)
La courbe ci dessus présente une saisonnalité annuelle avec une tendance stable
data_est_day %>% gg_season(Consommation, period = "year", labels = "both") +
labs(y="MW", title="Electricity demand: Grand Est")
data_est_day %>% gg_season(Consommation, period = "week", labels = "both") +
labs(y="MW", title="Electricity demand: Grand Est")
ggplot(data_est_day, aes(y=Consommation,x=TAVG))+
geom_point()+
theme_classic()+
geom_smooth(method="lm",
colour="red",
formula=y~x+I(x^2))
Nous voyons bien ici une corrélation polynomiale de type y = x + x² entre la température et la consommation électrique
## Visualisation
scatterplot(Consommation~dju_calc_total,data=data_est_day)
Idem pour ce graphique, nous notons une relation linéaire entre la consommation et les DJU
Ici, nous allons crée une variable consommation corrigé grâce aux régressions
lm_temp_corr<-lm(data_est_day$Consommation~data_est_day$dju_calc_total)
summary(lm_temp_corr)
##
## Call:
## lm(formula = data_est_day$Consommation ~ data_est_day$dju_calc_total)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1818.6 -406.6 102.1 427.4 1504.1
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4166.502 18.989 219.42 <2e-16 ***
## data_est_day$dju_calc_total 141.029 2.004 70.37 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 557.5 on 2554 degrees of freedom
## Multiple R-squared: 0.6597, Adjusted R-squared: 0.6596
## F-statistic: 4952 on 1 and 2554 DF, p-value: < 2.2e-16
a<-coef(lm_temp_corr)[2]
b<-coef(lm_temp_corr)[1]
Nous vérifions maintenant les conditions de validité de la régression linéaire
## Independance des résidus
acf(residuals(lm_temp_corr), main="lm_temp_corr")
durbinWatsonTest(lm_temp_corr)
## lag Autocorrelation D-W Statistic p-value
## 1 0.6127372 0.7718504 0
## Alternative hypothesis: rho != 0
# H0 (null hypothesis): There is no correlation among the residuals.
# HA (alternative hypothesis): The residuals are autocorrelated.
Nous voyons une forte autocorrelation des résidus sur le graphe et le test de Durbin Watson nous pousse à rejeter l’hypothèse de non corrélation
## Evaluation de l’hypothèse de normalité des résidus
plot(lm_temp_corr,2)
shapiro.test(residuals(lm_temp_corr))
##
## Shapiro-Wilk normality test
##
## data: residuals(lm_temp_corr)
## W = 0.97373, p-value < 2.2e-16
La normalité des résidus est accepté
## Evaluation de l’hypothèse d’homogénéité des résidus
plot(lm_temp_corr, 3)
ncvTest(lm_temp_corr)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 3.916416, Df = 1, p = 0.047817
L’homoscédasticité des residus est accepté
Les conditions de validité du modèle sont validées
data_est_day <-
data_est_day %>% mutate(conso_corrige_dju = Consommation - a * dju_calc_total)
data_est_month <-
data_est_month %>% mutate(conso_corrige_dju = Consommation - a * dju_calc_total)
library(xts)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following object is masked from 'package:tsibble':
##
## index
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
##
## Attaching package: 'xts'
## The following objects are masked from 'package:dplyr':
##
## first, last
data <- xts(x=data_est_day$conso_corrige_dju,order.by = data_est_day$Date_only)
dygraph(data)%>% dyRangeSelector()
data_est_month %>%
autoplot(conso_corrige_dju)
lm_temp_corr2<-lm(data_est_day$Consommation~data_est_day$TAVG+I(data_est_day$TAVG^2))
summary(lm_temp_corr2)
##
## Call:
## lm(formula = data_est_day$Consommation ~ data_est_day$TAVG +
## I(data_est_day$TAVG^2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1787.4 -414.3 110.0 387.3 1575.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6793.0822 24.7485 274.49 <2e-16 ***
## data_est_day$TAVG -189.4609 4.4057 -43.00 <2e-16 ***
## I(data_est_day$TAVG^2) 3.5738 0.1765 20.25 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 524.6 on 2553 degrees of freedom
## Multiple R-squared: 0.6989, Adjusted R-squared: 0.6986
## F-statistic: 2963 on 2 and 2553 DF, p-value: < 2.2e-16
x<-coef(lm_temp_corr2)[2]
x2<-coef(lm_temp_corr2)[3]
Nous vérifions maintenant les conditions de validité de la régression linéaire
## Indépendance des residus
acf(residuals(lm_temp_corr2), main="lm_temp_corr")
durbinWatsonTest(lm_temp_corr2)
## lag Autocorrelation D-W Statistic p-value
## 1 0.594772 0.8073433 0
## Alternative hypothesis: rho != 0
# H0 (null hypothesis): There is no correlation among the residuals.
# HA (alternative hypothesis): The residuals are autocorrelated.
Nous voyons une forte autocorrelation des résidus sur le graphe et le test de Durbin Watson nous pousse à rejeter l’hypothèse de non corr&lation
## Evaluation de l’hypothèse de normalité des résidus
plot(lm_temp_corr2,2)
shapiro.test(residuals(lm_temp_corr2))
##
## Shapiro-Wilk normality test
##
## data: residuals(lm_temp_corr2)
## W = 0.97537, p-value < 2.2e-16
La normalité des résidus est accepté
## Evaluation de l’hypothèse d’homogénéité des résidus
plot(lm_temp_corr2, 3)
ncvTest(lm_temp_corr2)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 16.43193, Df = 1, p = 5.0429e-05
L’homoscédasticité des résidus est accepté
Les conditions de validité du modèle sont validé
data_est_day <-
data_est_day %>% mutate(conso_corrige_temp = Consommation - abs((x * TAVG + (x2^2)*TAVG)))
data_temp <- data_est_day %>%
index_by(yearmonth(Date_only)) %>%
summarise(sum(conso_corrige_temp))
## merge des dataframmes
data_est_month <- left_join(data_est_month, data_temp, by = 'yearmonth(Date_only)')
library(xts)
data <- xts(x=data_est_day$conso_corrige_temp,order.by =data_est_day$Date_only)
dygraph(data)%>% dyRangeSelector()
Afin d’améliorer la compréhention et la précision du modèle, nous pouvons transformer celui ci.Soit en désaisonnalisant ou en appliquant une transformation de box-cox
Il s’agit d’une des premières méthode de desaisonnalisation, utilisé depuis 1920. Pour estimer la tendance mensuel avec une saisonnalité annuel, nous utilisons une moyenne mobile 2x12MA
data_est_month_ma <- data_est_month %>%
mutate(
`12-MA` = slider::slide_dbl(Consommation, mean,
.before = 5, .after = 6, .complete = TRUE),
`2x12-MA` = slider::slide_dbl(`12-MA`, mean,
.before = 1, .after = 0, .complete = TRUE)
)
data_est_month_ma %>%
autoplot(Consommation, colour = "gray") +
geom_line(aes(y = `2x12-MA`), colour = "#D55E00") +
labs(y = "MW",
title = "Consommation")
Une autre méthode utilisé par le census bureau est la méthode STL pour les séries additive sans variation calendaire, son principal avantage est de pouvoir traiter tout type de saisonnalité (heures, jours, mois, autres… )
## Desaisonnalisation de Consommation avec STL
dcmp_stl_conso <- data_est_month %>%
model(stl = STL(Consommation))
components(dcmp_stl_conso) %>%
as_tsibble() %>%
autoplot(Consommation, color="gray") +
geom_line(aes(y=trend+remainder), color = "#D55E00") +
labs(
y = "MV",
title = "electricity consommation"
)
components(dcmp_stl_conso) %>% autoplot()
## Desaisonnalisation de TAVG
dcmp_stl_tavg <- data_est_month %>%
model(stl = STL(TAVG))
components(dcmp_stl_tavg) %>%
as_tsibble() %>%
autoplot(TAVG, color="gray") +
geom_line(aes(y=trend), color = "#D55E00") +
labs(
y = "MV",
title = "Température moyenne journalière"
)
components(dcmp_stl_tavg) %>% autoplot()
## Desaisonnalisation de la conso corrigé avec STL
dcmp_conso <- data_est_day %>%
model(stl = STL(conso_corrige_dju))
components(dcmp_conso) %>%
as_tsibble() %>%
autoplot(conso_corrige_dju, color="gray") +
geom_line(aes(y=trend), color = "#D55E00") +
labs(
y = "MV",
title = "electricity consommation"
)
components(dcmp_conso) %>% autoplot()
## Ajout colonne adjusted au df
data_est_day$conso_corrige_adjusted <- components(dcmp_conso)$season_adjust
## Desaisonnalisation de TAVG
dcmp_tavg <- data_est_day %>%
model(stl = STL(TAVG))
components(dcmp_tavg) %>%
as_tsibble() %>%
autoplot(TAVG, color="gray") +
geom_line(aes(y=trend+remainder), color = "#D55E00") +
labs(
y = "MV",
title = "Température moyenne journalière"
)
components(dcmp_tavg) %>% autoplot()
Une méthode plus récente proposé par la Spain Bank est ARIMA SEATS X13, elle s’applique uniquement aux données mensuelles
dcmp_seats_conso <- data_est_month %>%
model(seats = X_13ARIMA_SEATS(conso_corrige_dju ~ seats())) %>%
components()
autoplot(dcmp_seats_conso) +
labs(title =
"Decomposition of conso_corrigé using SEATS")
dcmp_seats_conso %>%
ggplot(aes(x = `yearmonth(Date_only)`)) +
geom_line(aes(y = conso_corrige_dju, colour = "Data")) +
geom_line(aes(y = season_adjust,
colour = "Seasonally Adjusted")) +
geom_line(aes(y = trend, colour = "Trend")) +
labs(y = "MW",
title = "conso_corrigé (avg/month)") +
scale_colour_manual(
values = c("gray", "#0072B2", "#D55E00"),
breaks = c("Data", "Seasonally Adjusted", "Trend")
)
## Ajout colonne adjusted au df
data_est_month$conso_corrige_adjusted <- dcmp_seats_conso$season_adjust
dcmp_seats_tavg <- data_est_month %>%
model(seats = X_13ARIMA_SEATS(TAVG ~ seats())) %>%
components()
autoplot(dcmp_seats_tavg) +
labs(title =
"Decomposition of Temperature (avg/month) using SEATS")
dcmp_seats_tavg %>%
ggplot(aes(x = `yearmonth(Date_only)`)) +
geom_line(aes(y = TAVG, colour = "Data")) +
geom_line(aes(y = season_adjust,
colour = "Seasonally Adjusted")) +
geom_line(aes(y = trend, colour = "Trend")) +
labs(y = "°C ",
title = "Temperature (avg/month)") +
scale_colour_manual(
values = c("gray", "#0072B2", "#D55E00"),
breaks = c("Data", "Seasonally Adjusted", "Trend")
)
Nous percevons bien ici l’impact du réchauffement climatique, avec une montée de 2° en moyenne sur 7 ans
Nous pouvons utiliser une transformation de boxCox sur notre variable réponse afin normaliser les résidus. Nous choississons le parametre \(\\lambda\) grâce a l’indice de Guerrero
Cette transformation est utilisable directement dans fable, qui s’occupera alors de la transformation retour
## Transformation de Box Cox avec indice de Guerrero
lambda_guerrero_day <- data_est_month %>%
features(Consommation, features = guerrero) %>%
pull(lambda_guerrero)
data_est_month %>%
autoplot(box_cox(Consommation, lambda_guerrero_day)) +
#geom_line(aes(y=Consommation), colour = "#D55E00") +
labs(y = "",
title = latex2exp::TeX(paste0(
"Transformed electricity consommation with $\\lambda$ = ",
round(lambda_guerrero_day,2))))
gamma_guerrero_day <- data_est_month %>%
features(TAVG, features = guerrero) %>%
pull(lambda_guerrero)
data_est_month %>%
autoplot(box_cox(TAVG, gamma_guerrero_day)) +
geom_line(aes(y=TAVG), colour = "#D55E00") +
labs(y = "",
title = latex2exp::TeX(paste0(
"Transformed température with $\\lambda$ = ",
round(gamma_guerrero_day,2))))
# Estimate parameters
fit <- data_est_month %>%
model(ETS(Consommation ~ error("A") + trend("N") + season("N")))
fc <- fit %>%
forecast(h ="1 years")
fc %>%
autoplot(data_est_month) +
geom_line(aes(y = .fitted), col="#D55E00",
data = augment(fit)) +
labs(y="MW", title="Consommation") +
guides(colour = FALSE)
# Estimate parameters
fit <- data_est_month %>%
model(ETS(Consommation ~ error("A") + trend("A") + season("N")))
fc <- fit %>%
forecast(h ="1 years")
fc %>%
autoplot(data_est_month) +
geom_line(aes(y = .fitted), col="#D55E00",
data = augment(fit)) +
labs(y="MW", title="Consommation") +
guides(colour = FALSE)
# Estimate parameters
fit <- data_est_month %>%
model(ETS(Consommation ~ error("A") + trend("A") + season("A")))
fc <- fit %>%
forecast(h ="1 years")
fc %>%
autoplot(data_est_month) +
geom_line(aes(y = .fitted), col="#D55E00",
data = augment(fit)) +
labs(y="MW", title="Consommation") +
guides(colour = FALSE)
Testons enfin un dernier modèle en full auto
# Estimate parameters
fit <- data_est_month %>%
model(ETS(Consommation))
fc <- fit %>%
forecast(h ="1 years")
report(fit)
## Series: Consommation
## Model: ETS(M,N,M)
## Smoothing parameters:
## alpha = 0.003656971
## gamma = 0.0002115003
##
## Initial states:
## l s1 s2 s3 s4 s5 s6 s7
## 160457.7 1.17136 1.098812 0.9908194 0.8667747 0.7817055 0.8653705 0.8523303
## s8 s9 s10 s11 s12
## 0.9021363 0.9519064 1.12806 1.13257 1.258155
##
## sigma^2: 0.0013
##
## AIC AICc BIC
## 1839.675 1846.734 1876.137
fc %>%
autoplot(data_est_month) +
geom_line(aes(y = .fitted), col="#D55E00",
data = augment(fit)) +
labs(y="MW", title="Consommation") +
guides(colour = FALSE)
AR = Auto regressive
I = Integrated
MA = Mooving average
S = Seaseonal
Approche plutôt basé sur l’autocorrélation
Nous n’utilisons plus les Trend et saisonnalité pour modéliser le processus, nous devons donc stationnarisé la série temporelle, la rendre independante du temps.
Nous regardons dans un premier temps si nous avons besoin de différencier grâce a un test de racine unitaire
data_est_month %>%
features(Consommation,unitroot_nsdiffs)
## # A tibble: 1 x 1
## nsdiffs
## <int>
## 1 1
Le test saisonnier nous préconise de différencier 1 fois la variable consommation en saisonnalité
Faisons le même test mais avec le test en tendance
data_est_month %>%
features(Consommation, unitroot_ndiffs)
## # A tibble: 1 x 1
## ndiffs
## <int>
## 1 0
Analysons le graphes des résidus de la consommation différencié
data_est_month %>%
gg_tsdisplay(difference((Consommation),12), plot_type='partial')
Recommençons pour la variable conso désaisonnalisé.
data_est_month %>%
features(conso_corrige_adjusted,unitroot_nsdiffs)
## # A tibble: 1 x 1
## nsdiffs
## <int>
## 1 0
data_est_month %>%
features(conso_corrige_adjusted,unitroot_ndiffs)
## # A tibble: 1 x 1
## ndiffs
## <int>
## 1 0
Le test confirme que cette variables est déjà differencié en tendance et saisonnalité
data_est_month %>%
gg_tsdisplay((conso_corrige_adjusted), plot_type='partial')
Nous utilisons la variable Consommation corrigé de l’effet température et désaisonnalisé
fit_arima <- data_est_month %>%
model(
auto = ARIMA(conso_corrige_adjusted),
arima012011 = ARIMA(conso_corrige_adjusted ~ pdq(3, 0, 2) + PDQ(0, 0, 0))## PDQ a 0 pour indiquer qu'il n'y a pas de tendance saisonière
# model_search = ARIMA(conso_corrige_dju ~ pdq(p=0:2, d=1, q=0:2) + PDQ(0,1,1))
)
forecast(fit_arima,h=12) %>%
autoplot(data_est_month)
fit_arima %>% pivot_longer(everything(), names_to = "Model name",values_to = "Orders")
## # A mable: 2 x 2
## # Key: Model name [2]
## `Model name` Orders
## <chr> <model>
## 1 auto <ARIMA(1,0,1) w/ mean>
## 2 arima012011 <ARIMA(3,0,2) w/ mean>
glance(fit_arima) %>% arrange(AICc) %>% select(.model:BIC)
## # A tibble: 2 x 6
## .model sigma2 log_lik AIC AICc BIC
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 auto 6303449. -775. 1559. 1559. 1568.
## 2 arima012011 6407252. -774. 1563. 1564. 1580.
fit_arima %>% select(auto) %>% report()
## Series: conso_corrige_adjusted
## Model: ARIMA(1,0,1) w/ mean
##
## Coefficients:
## ar1 ma1 constant
## 0.7702 -0.5465 29134.9501
## s.e. 0.1291 0.1586 119.0057
##
## sigma^2 estimated as 6303449: log likelihood=-775.32
## AIC=1558.64 AICc=1559.15 BIC=1568.37
On remarque que le modèle ne performe pas très bien du fait que notre variable ressemble de plus en plus a du bruit blanc et donc très difficile de prédire
| p=p= | order of the autoregressive part; |
|---|---|
| d=d= | degree of first differencing involved; |
| q=q= | order of the moving average part. |
fit_sarima <- data_est_month %>%
model(
auto = ARIMA(conso_corrige_dju),
arima012011 = ARIMA(conso_corrige_dju ~ pdq(0, 0, 2) + PDQ(0, 1, 1))
# model_search = ARIMA(conso_corrige_dju ~ pdq(p=0:2, d=1, q=0:2) + PDQ(0,1,1))
)
forecast(fit_sarima,h=24) %>%
autoplot(data_est_month)
fit_sarima %>% pivot_longer(everything(), names_to = "Model name",values_to = "Orders")
## # A mable: 2 x 2
## # Key: Model name [2]
## `Model name` Orders
## <chr> <model>
## 1 auto <ARIMA(1,0,1)(2,1,0)[12]>
## 2 arima012011 <NULL model>
glance(fit_sarima) %>% arrange(AICc) %>% select(.model:BIC)
## # A tibble: 1 x 6
## .model sigma2 log_lik AIC AICc BIC
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 auto 12303309. -691. 1391. 1392. 1402.
Nous devons maintenant réintroduire la correction de température
Nous pouvons maintenant comparé les modèle grace a la cross-validation
data_est_month_tr <- data_est_month %>%
stretch_tsibble(.init = 3, .step = 1)
fc <- data_est_month_tr %>%
model(
auto = ARIMA(conso_corrige_dju)
) %>%
forecast(h = 24) %>%
group_by(.id) %>%
mutate(h = row_number()) %>%
ungroup()
fc %>%
accuracy(data_est_month, by = c("h", ".model")) %>%
ggplot(aes(x = h, y = MAPE)) +
geom_point()
Ici, nous ajoutons directement des variables éxogène à notre modèle
training <- data_est_month %>% filter(year(`yearmonth(Date_only)`) <= 2018)
test <- data_est_month %>% filter(year(`yearmonth(Date_only)`) > 2018)
fit_sarimax <- data_est_month %>%
model(
auto = ARIMA(Consommation~ TAVG + I(TAVG^2)), ## Auto
sarimax101210 = ARIMA(Consommation~ TAVG + I(TAVG^2)+pdq(1,0,1)+PDQ(2,1,0)) ## idem modele correction
)
fit_sarimax %>% pivot_longer(everything(), names_to = "Model name",values_to = "Orders")
## # A mable: 2 x 2
## # Key: Model name [2]
## `Model name` Orders
## <chr> <model>
## 1 auto <LM w/ ARIMA(0,0,1)(2,1,0)[12] errors>
## 2 sarimax101210 <LM w/ ARIMA(1,0,1)(2,1,0)[12] errors>
glance(fit_sarimax) %>% arrange(AICc) %>% select(.model:BIC)
## # A tibble: 2 x 6
## .model sigma2 log_lik AIC AICc BIC
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 auto 10431952. -685. 1382. 1383. 1396.
## 2 sarimax101210 10538566. -685. 1383. 1385. 1399.
fit_sarimax %>%
select(auto) %>%
report()
## Series: Consommation
## Model: LM w/ ARIMA(0,0,1)(2,1,0)[12] errors
##
## Coefficients:
## ma1 sar1 sar2 TAVG I(TAVG^2)
## 0.4284 -0.6336 -0.4129 -4795.0171 119.6446
## s.e. 0.1292 0.1225 0.1234 317.9231 14.0748
##
## sigma^2 estimated as 10431952: log likelihood=-685.01
## AIC=1382.03 AICc=1383.32 BIC=1395.69
fit_sarimax %>% accuracy()
## # A tibble: 2 x 10
## .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 auto Training -360. 2885. 2087. -0.254 1.31 0.297 0.286 -0.0418
## 2 sarimax101210 Training -386. 2878. 2105. -0.272 1.33 0.300 0.285 -0.00312
# Time series cross-validation accuracy
data_est_month_tr <- data_est_month %>%
stretch_tsibble(.init = 1, .step = 12)
data_est_month_tr
## # A tsibble: 259 x 14 [1M]
## # Key: .id [7]
## `yearmonth(Date_only)` TAVG TMAX TMIN Consommation Renouvelable
## <mth> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2013 janv. 2.25 4.33 0.103 202862. 42650.
## 2 2013 janv. 2.25 4.33 0.103 202862. 42650.
## 3 2013 févr. 0.968 3.71 -1.45 191247. 41756.
## 4 2013 mars 3.56 7.56 0.103 195760. 40565.
## 5 2013 avr. 10.9 16.0 6.30 161972. 48322.
## 6 2013 mai 12.9 17.2 9.05 150792. 48971.
## 7 2013 juin 18.7 24.0 13.2 137521. 48742.
## 8 2013 juil. 22.4 28.3 15.9 140700. 46800.
## 9 2013 août 19.7 26.2 13.8 124791. 34938
## 10 2013 sept. 15.8 21.5 11.3 141255. 37687.
## # ... with 249 more rows, and 8 more variables: dju_calc_chauffage <dbl>,
## # dju_calc_clim <dbl>, dju_calc_total <dbl>, PRCP <dbl>,
## # conso_corrige_dju <dbl>, sum(conso_corrige_temp) <dbl>,
## # conso_corrige_adjusted <dbl>, .id <int>
# TSCV accuracy
fit <- data_est_month_tr %>%
model(ARIMA(Consommation)) %>%
forecast(h = "2 years") %>%
accuracy(data_est_month)
forecast(fit_sarimax,new_data = test) %>%
autoplot(data_est_month) +
labs(title = "Electricity demand",
y = "MW") + coord_cartesian(xlim = c(as_date("2018-01-01"),as_date("2019-12-31")))
Cette méthode est basé sur les réseaux de neurones artificiels, nous pouvons aussi y ajouter des variables explicatives.
Le temps de calcul étant relativement long, nous travaillons sur le data month. Vu que nous utilisons une variable exogene dans notre modèle, nous splittons le dataset en Train/Test
fit_mnn <- training %>%
model(
AutoX = NNETAR(Consommation , xreg = training$TAVG),
Auto = NNETAR(Consommation)
)
forecast(fit_mnn,new_data = test, times = 100) %>% ## times = 200 permets de limiter le temps de calcul en dimuniant la précisions des intervalles de confiance
autoplot(data_est_month) +
labs(title = "Electricity demand",
y = "MW") + coord_cartesian(xlim = c(as_date("2018-01-01"),as_date("2019-12-31"))) ## Restiction pour ne voir que 2 ans de 2018 a 2019
accuracy(fit_mnn)
## # A tibble: 2 x 10
## .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 AutoX Training -4.62 7281. 5219. -0.191 3.09 0.747 0.715 0.348
## 2 Auto Training 8.03 7366. 5237. -0.186 3.10 0.750 0.724 0.334
Les résultats sont convainquant avec un MAPE de 3%
Le modèle Prohet est un modèle récent de 2018 (Taylor and Letham 2018) viens directement de Facebook, basé sur des régressions avec des séries de Fourier.
Ses avantage sont sa simplicité d’utilisation et son intégration des effets vacances, l’incovenient principale est que nous ne maitrisons pas les choix de modèle, basé sur des critère Bayésien
library(prophet)
## Loading required package: Rcpp
## Loading required package: rlang
##
## Attaching package: 'rlang'
## The following objects are masked from 'package:purrr':
##
## %@%, as_function, flatten, flatten_chr, flatten_dbl, flatten_int,
## flatten_lgl, flatten_raw, invoke, list_along, modify, prepend,
## splice
library(dygraphs)
Nous importons un fichiers avec les jours fériés Francais, nous pourrions importer les dates des vacances mais cela nécessite plus de manipulation. De plus, ils devrais être integrés prochainement directement dans Prophet
jours_feries <- read.csv("data/holidays_fra.csv",fileEncoding="UTF-8")
Nous devons renomer les variables date (ds) et d’interet (y)
data_prophet_day <- data_est_day %>%
rename(ds = Date_only, y = Consommation)
data_prophet_month <- data_est_month %>%
rename(ds = `yearmonth(Date_only)`, y = Consommation)
Nous renseignons le modele
m<- prophet(holidays = jours_feries) ## ajout des jours feries
m<- add_regressor(m,'TAVG') ## ajout variables expliquatives
m<- fit.prophet(m,data_prophet_day) ## fit du model
## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.
Nous devons crées un dataframme future a partir duquel nous pourrons faire les prévision. Pour le moment, nous avons que les dates dedans, il nous faut les prévisions de la variable exliquative. Nous pourrions l’importer depuis un organisme météo, mais j’ai choisi de la prévoir directement à l’aide de Prophet
data_prophet_temp2 <- data_prophet_day %>%
rename(Consomation=y,y = TAVG) ## renommage des variables
m2 <- prophet(data_prophet_temp2) ## on renseigne le df
## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.
future <- make_future_dataframe(m2, periods = 365) ## on crée un df future de 1 an
forecast <- predict(m2, future) ## on realise les predictions
forecast_temp_2020 <- data.frame(forecast[c('ds', 'yhat')]) %>% ## on transforme les prédictions en df
rename(TAVG=yhat) ## on renome TAVG
Nous pouvons maintenant lancer les prévisions.
forecast <- predict(m, forecast_temp_2020) ## predictions
prophet_plot_components(m, forecast)## composante
dyplot.prophet(m, forecast,uncertainty = TRUE)## Affichage dynamique
Vérifions maintenant la qualité des prédictions à l’aide d’une cross validation
df.cv <- cross_validation(m, horizon=365, units='days')
## Making 7 forecasts with cutoffs between 2016-01-01 and 2018-12-31
plot_cross_validation_metric(df.cv, metric='mape')
Nous obtenons un résultat très intéressant, avec une erreur moyenne de 5%